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Abstract 

We provide a complete picture of the non-relativistic self-gravitating gas at thermal equilibrium 
using Monte Carlo simulations, analytic mean field methods (MF) and low density expansions. The 
system is shown to possess an infinite volume limit in the grand canonical (GCE), canonical (CE) 
and microcanonical (MCE) ensembles when (N, V) — > oo, keeping N/V 1 ^ 3 fixed. We compute the 
equation of state (we do not assume it as is customary), as well as the energy, free energy, entropy, 
chemical potential, specific heats, compressibilities and speed of sound; we analyze their properties, 
signs and singularities. All physical quantities turn out to depend on a single variable r\ = y™ 3 ^ 
that is kept fixed in the N — > oo and V —* oo limit. The system is in a gaseous phase for r\ <t\t and 
collapses into a dense object for rj > t]t in the CE with the pressure becoming large and negative. 
At rj = t)t the isothermal compressibility diverges; this gravitational phase transition is associated 
to the Jeans' instability, our Monte Carlo simulations yield rjr — 1.515. PV/[NT] — f(rf) and all 
physical magnitudes exhibit a square root branch point at rj = rjc > t]t. The values of ryr and r/c 
change by a few percent with the geometry for large N: for spherical symmetry and N = oo (MF), 
we find rjc = 1.561764 . . ., while the Monte Carlo simulations for cubic geometry yields rjc — 1.540. 
In mean field and spherical symmetry cv diverges as (rjc — ?7)~ 1/ ' 2 f° r V T Vc while cp and kt diverge 
as (t)t — v)^ 1 f° r V T Vt — 1.51024 . . .. The function f(rj) has a second Riemann sheet which is only 
physically realized in the MCE. In the MCE, the collapse phase transition takes place in this second 
sheet near tjmc = 1-26 and the pressure and temperature are larger in the collapsed phase than in 
the gaseous phase. Both collapse phase transitions (in the CE and in the MCE) are of zeroth order 
since the Gibbs free energy has a jump at the transitions. The MF equation of state in a sphere, 
f(rj), obeys a first order non-linear differential equation of first kind Abel's type. The MF gives 
an extremely accurate picture in agreement with the MC simulations both in the CE and MCE. We 
perform the MC simulations on a cubic geometry, which thus describe an isothermal cube while 
the MF calculations describe an isothermal sphere. 

We complete our study of the self-gravitating gas by computing the fluctuations around the sad- 
dle point solution for the three statistical ensembles (grand canonical, canonical and microcanonical). 
Although the saddle point is the same for the three ensembles, the fluctuations change from one 
ensemble to the other. The zeroes of the small fluctuations determinant determine the position of 
the critical points for each ensemble. This yields the domains of validity of the mean field approach. 
Only the S wave determinant exhibits critical points. Closed formulae for the S and P wave determi- 
nants of fluctuations are derived. The local properties of the self-gravitating gas in thermodynamic 
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equilibrium are studied in detail. The pressure, energy density, particle density and speed of sound 
are computed and analyzed as functions of the position. The equation of state turns out to be locally 
p(r) = T pv(f) as for the ideal gas. Starting from the partition function of the self-gravitating gas, 
we prove in this microscopic calculation that the hydrostatic description yielding locally the ideal 
gas equation of state is exact in the N = oo limit. The dilute nature of the thermodynamic limit 
(N ~ L — > oo with N/L fixed) together with the long range nature of the gravitational forces play a 
crucial role in obtaining such ideal gas equation. The self-gravitating gas being inhomogeneous, we 
have PV/[NT] = f(rj) < 1 for any finite volume V. The inhomogeneous particle distribution in the 
ground state suggests a fractal distribution with Haussdorf dimension D, D is slowly decreasing with 
increasing density, 1 < D < 3. The average distance between particles is computed in Monte Carlo 
simulations and analytically in the mean field approach. A dramatic drop at the phase transition is 
exhibited, clearly illustrating the properties of the collapse. 
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1 Statistical Mechanics of the Self- Gravitating Gas 

Physical systems at thermal equilibrium are usually homogeneous. This is the case for gases with short 
range intermolecular forces (and in absence of external fields). In such cases the entropy is maximum 
when the system homogenizes. 

When long range interactions as the gravitational force are present, even the ground state is inho- 
mogeneous. In this case, each element of the substance is acted on by very strong forces due to distant 
particles in the gas. Hence, regions near to and far from the boundary of the volume occupied by the gas 
will be in very different conditions, and, as a result, the homogeneity of the gas is destroyed py. The state 
of maximal entropy for gravitational systems is inhomogeneous. This basic inhomogeneity suggested 
us that fractal structures can arise in a self- interacting gravitational gas 0121 EB • 

The inhomogeneous character of the ground state for gravitational systems explains why the universe 
is not going towards a 'thermal death'. A 'thermal death' would mean that the universe evolves towards 
more and more homogeneity. This can only happen if the entropy is maximal for an homogeneous state. 
Instead, it is the opposite what happens, structures are formed in the universe through the action of the 
gravitational forces as time evolves. 

Usual theorems in statistical mechanics break down for inhomogeneous ground states. For example, 
the specific heat may be negative in the microcanonical ensemble (not in the canonical ensemble where 
it is always positive) 3 . 

As is known, the thermodynamic limit for self-gravitating systems does not exist in its usual form 
(N — > 00, V — » 00, N/V = fixed). The system collapses into a very dense phase which is determined 
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by the short distance (non-gravitational) forces between the particles. However, the thermodynamic 
functions exist in the dilute limit ^ |2] 



N -> oo , V 



N 



— fixed , 



oo , 



where V stands for the volume of the box containing the gas. In such a limit, the energy E, the free 
energy and the entropy turns to be extensive. That is, we find that they take the form of N times a 
function of the intensive dimcnsionless variables: 



where rj and £ are intensive variables. Namely, r\ and £ stay finite when N and V = L 3 tend to infinite. 
The variable rj is appropriate for the canonical ensemble and £ for the microcanonical ensemble. Physical 
magnitudes as the specific heat, speed of sound, chemical potential and compressibility only depend on 
77 or £. The variables r\ and £, as well as the ratio N/L, are therefore intensive magnitudes. The energy, 
the free energy, the Gibbs free energy and the entropy are of the form TV times a function of r\. These 
functions of 77 have a finite N — 00 limit for fixed rj (once the ideal gas contributions are subtracted). 
Moreover, the dependence on rj in all these magnitudes express through a single universal function f{rf). 
The variable r\ is the ratio of the characteristic gravitational energy G " — and the kinetic energy T of a 
particle in the gas. For 77 = the ideal gas is recovered. 

In refs. ^G] we have thoroughly studied the statistical mechanics of the self-gravitating gas. That 
is, our starting point is the partition function for non-relativistic particles interacting through their 
gravitational attraction in thermal equilibrium. We study the self-gravitating gas in the three ensembles: 
microcanonical (MCE), canonical (CE) and grand canonical (GCE). We performed calculations by three 
methods: 

• By expanding the partition function through direct calculation in powers of 1 /£ and rj for the MCE 
and CE, respectively. These expressions apply in the dilute regime (£ 3> 1 , rj <C 1) and become 
identical for both ensembles for N — > 00. At 77 = 01/£ we recover the ideal gas behaviour. 

• By performing Monte Carlo simulations both in the MCE and in the CE. We found in this way that 
the self-gravitating gas collapses at a critical point which depends on the ensemble considered. As 
shown in fig. Q]the collapse occurs first in the canonical ensemble (point T). The microcanonical 
ensemble exhibits a larger region of stability that ends at the point MC (fig. Notice that the 
physical magnitudes are identical in the common region of validity of both ensembles within the 
statistical error. Beyond the critical point T the system becomes suddenly extremely compact 
with a large negative pressure in the CE. Beyond the point MC in the MCE the pressure and the 
temperature increase suddenly and the gas collapses. The phase transitions at T and at MC are 
of zeroth order since the Gibbs free energy has discontinuities in both cases. 

• By using the mean field approach we evaluate the partition function for large N . We do this 
computation in the grand canonical, canonical and microcanonical ensembles. In the three cases, 
the partition function is expressed as a functional integral over a statistical weight which depends 
on the (continuous) particle density. These statistical weights are of the form of the exponential of 
an 'effective action' proportional to N . Therefore, the N — > 00 limit follows by the saddle point 
method. The saddle point is a space dependent mean field showing the inhomogeneous character 
of the ground state. Corrections to the mean field are of the order 1/N and can be safely ignored 
for N 3> 1 except near the critical points. These mean field results turned out to be in excellent 
agreement with the Monte Carlo results and with the low density expansion. 

We calculate the saddle point (mean field) for spherical symmetry and we obtain from it the various 
physical magnitudes (pressure, energy, entropy, free energy, specific heats, compressibilities, speed of 
sound and particle density). Furthermore, we computed in ref.[2] the determinants of small fluctua- 
tions around the saddle point solution for spherical symmetry for the three statistical ensembles. 

When any small fluctuation around the saddle point decreases the statistical weight in the functional 
integral, the saddle point dominates the functional integral and the mean field approach can be valid. In 
that case, the determinant of small fluctuations is positive. A negative determinant of small fluctuations 
indicates that some fluctuations around the saddle point are increasing the statistical weight in the 
functional integral and hence the saddle point does not dominate the partition function. The mean field 
approach cannot be used when the determinant of small fluctuations is negative. We find analytically 
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Figure 1: f(rj R ) = PV/[NT] as a function of r] R by Monte Carlo simulations for the microcanonical and 
canonical ensembles (N = 2000). Both curves coincide within the statistical error till the point T. 

in the CE that the determinant of small fluctuations vanishes at the point ijc = 1.561764 and becames 
negative at 77 > r\c Qj- (The point r\c is indicated C in figCJ. We find that the CE specific heat cy at 
the point C diverges as PQ: 



The (+) sign refers to the positive (first) branch, and the (-) sign to the negative (second) branch (between 
the points C and MC). 

However, it must be noticed that the instability point is located at rj = r)x < r)c, as shown by 
both, mean field and Monte Carlo computations. (The point t\t is indicated T in fig. The onset 
of instability in the canonical ensemble coincides with the point where the isothermal compressibility 
diverges. The isothermal compressibility kt is positive from 77 = till 77 = Tfr = 1.51024. . .. At this 
point kt as well as the specific heat at constant pressure cp diverge and change their signs 1. Moreover, 
at this point the speed of sound at the center of the sphere becomes imaginary 2 . Therefore, small 
density fluctuations will grow exponentially in time instead of exhibiting oscillatory propagation. Such 
a behaviour leads to the collapse of the gas into a extremely compact object. Monte-Carlo simulations 
confirm the presence of this instability at t\t = 1-510 ... in the canonical ensemble and the formation of 
the collapsed object pQ. 

The collapse in the Grand Canonical ensemble (GCE) occurs for a smaller value of r\ — j]gc = 
0.49465 . . ., while in the microcanonical ensemble, the collapse arrives later, in the second sheet, at 
V = Vmc = 1.25984.... 

We find that the Monte Carlo simulations for self-gravitating gas in the CE and the MCE confirm 
the stability results obtained from mean field. 

The saddle point solution is identical for the three statistical ensembles. This is not the case for the 
fluctuations around it. The presence of constraints in the CE (on the number of particles) and in the 
MCE (on the energy and the number of particles) changes the functional integral over the quadratic 
fluctuations with respect to the GCE. 

The saddle point of the partition function turns out to coincide with the hydrostatic treatment of 
the self-gravitating gas ^3 E| (which is usually known as the 'isothermal sphere' in the spherically 
symmetric case). 

We find that the Monte Carlo simulations (describing thermal equilibrium) are much more efficient 
than the TV-body simulations integrating Newton's equations of motion. (Indeed, the integration of 



c v = ±0.63572 . . . (77c - V)~ 1/2 - 0.19924 . . . + 0(Vvc - v) 
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Microcanonical, collapsed, zi=-0.6 




Figure 2: Average particle distribution in the collapsed phase from Monte Carlo simulations with 2000 
particles in the microcanonical ensemble for £ = —0.6, T) = 0.43, pV/[NT] = 0.414. 



Newton's equations provides much more detailed information than the one needed in thermal equilibrium 
investigations). Actually, a few hundreds of particles are enough to get quite accurate results in the Monte 
Carlo simulations (except near the collapse points). Moreover, the Monte Carlo results turns to be in 
excellent agreement with the mean field calculations up to very small corrections of the order 1/N. Our 
Monte Carlo simulations are performed in a cubic geometry. The equilibrium configurations obtained in 
this manner can thus be called the 'isothermal cube'. 

In summary, the picture we get from our calculations using these three methods show that the self- 
gravitating gas behaves as a perfect gas for r\ — > 0, l/£ — > 0. When r\ and l/£ grow, the gas becomes 
denser till it suddenly condenses into a high density object at a critical point GC, C or MC depending 
upon the statistical ensemble chosen. 

r\ is related with the Jeans' length dj of the gas through r\ = 3 (L/dj) 2 . Hence, when rj goes beyond 
t]t, the length of the system becomes larger than dj / 'y 't\t /3. The collapse at T in the CE is therefore a 
manifestation of the Jeans' instability. 

In the MCE, the determinant of fluctuations vanishes at the point MC. The physical states beyond 
MC are collapsed configurations as shown by the Monte Carlo simulations [see fig. 0. Actually, the gas 
collapses in the Monte Carlo simulations slightly before the mean field prediction for the point MC. The 
phase transition at the microcanonical critical point MC is the so called gravothermal catastrophe |13| . 

The gravitational interaction being attractive without lower bound, a short distance cut-off (A) must 
be introduced in order to give a meaning to the partition function. We take the gravitational force 
between particles as — G m 2 /r 2 for r > A and zero for r < A where r is the distance between the two 
particles. We show that the cut-off effects are negligible in the N = oo limit. That is, once we set N = oo 
with fixed 77, all physical quantities are finite in the zero cut-off limit (A — 0). The cut-off effects are of 
the order A 2 /L 2 and can be safely ignored. 

In ref. |T| we expressed all global physical quantities in terms of a single function /(??). Besides 
computing numerically f(rf) in the mean field approach, we showed that this function obeys a first order 
non-linear differential equation of first Abel's type. We obtained analytic results about f(rj) from the 
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Abel's equation. f(rj) exhibits a square-root cut at rjc, the critical point in the CE. The first Riemann 
sheet is realized both in the CE and the MCE, whereas the second Riemann sheet (where cy < 0) is only 
realized in the MCE. f(rf) has infinitely many branches in the 77 plane but only the first two branches are 
physically realized. Beyond MC the states described by the mean field saddle point are unstable. 

We plot and analyze the equation of state, the energy, the entropy, the free energy, cy and the 
isothermal compressibility [figs. 161 lUj . Most of these physical magnitudes were not previously computed 
in the literature as functions of 77. 

We find analytically the behaviour of f(rj) near the point r\c i n mean field, 

(r)c-vf /2 ~ ■ 

This exhibits a square root branch point singularity at 77c- This shows that the specific heat at constant 

1/2 

volume diverges at rjc as (i]c — rj) for 77 f rjc- However, it must be noticed, that the specific heat at 
constant pressure and the isothermal compressibility both diverge at the point r/T as (t]t — rj) . These 
mean field results apply for \t] — tjc\ 1 <C N \ rj — rjc |- Fluctuations around mean field can be neglected 
in such a regime. 

The Monte Carlo calculations permit us to obtain f{ri) in the collapsed phase. Such result cannot be 
obtained in the mean field approach. The mean field approach only provides information as f(rj) in the 
dilute gas phase. 

For the self-gravitating gas, we find that the Gibbs free energy $ is not equal to N times the chemical 
potential and that the thermodynamic potential fl is not equal to — PV as usual [jj]. This is a consequence 
of the dilute thermodynamic limit N — > 00, L — > 00, N/L =fixed. 

We computed the determinant of small fluctuations around the saddle point solution for spherical 
symmetry in all three statistical ensembles. In the spherically symmetric case, the determinant of small 
fluctuations is written as an infinite product over partial waves. The S and P wave determinants are 
written in closed form in terms of the saddle solution. The determinants for higher partial waves are 
computed numerically. All partial wave determinants are positive definite except for the S-wave P_ . The 
reason why the fluctuations are different in the three ensembles is rather simple. The more contraints 
are imposed the smaller becomes the space of fluctuations. Therefore, in the grand canonical ensemble 
(GCE) the system is more free to fluctuate and the phase transition takes place earlier than in the micro- 
canonical (MCE) and canonical ensembles (CE). For the same reason, the transition takes place earlier 
in the CE than in the MCE. 

The conclusion being that the mean field correctly gives an excellent description of the thermodynamic 
limit except near the critical points (where the small fluctuations determinant vanishes); the mean field 
is valid for N\rj — r) cr u\ 1. The vicinity of the critical point should be studied in a double scaling limit 
N — > 00, rj — > rjcrit- Critical exponents are reported in ref. [lj for r\ — > rjc using the mean field. These 
mean field results apply for \rj — rjc\ ^ 1 ^ N\r) — rjc\ with N ^> 1. Fluctuations around mean field can 
be neglected in such a regime. 

We computed local properties of the gas in ref.fSj. That is, the local energy density e(r), local particle 
density, local pressure and the local speed of sound. Furthermore, we analyze the scaling behaviour of 
the particle distribution and its fractal (Haussdorf) dimension 

The particle distribution pv{<t) proves to be inhomogeneous (except for 77 -C 1) and described by 
an universal function of 77, the geometry and the ratio r = q/ R 1 R being the radial size. Both Monte 
Carlo simulations and the Mean Field approach show that the system is inhomogeneous forming a clump 
of size smaller than the box of volume V [see figs. I3I5| . 

The particle density in the bulk behaves as pv{<t) — r D ~ 3 . That is, the mass M{R) enclosed on a 
region of size R vary approximately as 

M(R) ~ C R D . 

D slowly decreases from the value D — 3 for the ideal gas (77 = 0) till D = 0.98 in the extreme limit of the 
MC point, D takes the value 1.6 at 77c, [see Table 2]. This indicates the presence of a fractal distribution 
with Haussdorf dimension D. 

Our study of the statistical mechanics of a self-gravitating system indicates that gravity provides a 
dynamical mechanism to produce fractal structures E] ■ 

The average distance between particles monotonically decrease with 77 in the first sheet. The mean 
field and Monte Carlo are very close in the gaseous phase whereas the Monte Carlo simulations exhibit a 
spectacular drop in the average particle distance at the clumping transition point T. In the second sheet 
(only described by the MCE) the average particle distance increases with 77 



f(v) y o + °- 27137 • • ■ Vv^V + 0.27763 . . . ( Vc - rj) + O 
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Wc find that the local equation of state is given by 



p{r) = Tp v (r) . 



(1) 



We have derived the equation of state for the self-gravitating gas. It is locally the ideal gas equation, 
but the self-gravitating gas being inhomogeneous, the pressure at the surface of a given volume is not 
equal to the temperature times the average density of particles in the volume. In particular, for the whole 
volume: PV/[NT] = f(rj) < 1 (the equality holds only for 77 = 0). 

Notice that we have found the local ideal gas equation of state p(r) — Tpy (r ) for purely gravitational 
interaction between particles. Therefore, equations of state different from this one, (as often assumed 
and used in the literature for the self-gravitating gas), necessarily imply the presence of additional non- 
gravitational forces. 

The local energy density e(r) turns out to be an increasing function of r in the spherically symmetric 
case. The energy density is always positive on the surface, whereas it is positive at the center for 
< rj < 773 = 1.07783 . . ., and negative beyond the point 77 = 773 = 1.07783 . . .. 

The local speed of sound v 2 s (r) is computed in the mean field approach as a function of the position 
for spherical symmetry and long wavelengths. Vg(r) diverges at 77 = rjr = 1.51024 ... in the first Riemann 
sheet. Just beyond this point uf(r) is large and negative in the bulk showing the strongly unstable 
behaviour of the gas for such range of values of 77. 

Moreover, we have shown the equivalence between the statistical mechanical treatment in the mean 
field approach and the hydrostatic description of the self-gravitating gas ^3 EI ■ 

The success of the hydrodynamical description depends on the value of the mean free path (I) compared 
with the relevant sizes in the system. I must be < 1. We compute the ratio I /a (Knudsen number), 
where a is a length scale that stays fixed for N — > 00 and show that l/a ~ iV~ 2 . This result ensures the 
accuracy of the hydrodynamical description for large N . 

Furthermore, we have computed in ref. |2] several physical magnitudes as functions of 77 and r 
which were not previously computed in the literature as the speed of sound, the energy density, the 
average distance between particles and we notice the presence of a Haussdorf dimension in the particle 
distribution. 

The statistical mechanics of a selfgravitating gas formed by particles with different masses is thor- 
oughly investigated in ref. [7] while the selfgravitating gas in the presence of the cosmological constant 
is thoroughly investigated in refs. [201 ■ In ref.[27j the Mayer expansion for the selfgravitating gas is 
investigated in connection with the stability of the gaseous phase. 

The variable rj R appropriate for a spherical symmetry is defined as 



We use 77 for the mean field theory calculations with spherical symmetry while 77 is preferred for the 
Monte Carlo simulations with cubic symmetry. 

2 Statistical Mechanics of the Self-Gravitating Gas: the micro- 
canonical and the canonical ensembles 

We investigate a gas of N non-relativistic particles with mass m self-interacting through Newtonian 
gravity. We consider first the particles isolated, that is, a self-gravitating gas in the microcanonical 
ensemble. We subsequently consider them in thermal equilibrium at temperature T = /3 _1 . That is, in 
the canonical ensemble where the system of N particles is not isolated but in contact with a thermal 
bath at temperature T. 

We always assume the system being on a cubic box of side L just for simplicity. We consider spherical 
symmetry for the mean field approach (sec. 4). Please notice that we never use periodic boundary 
conditions. 

At short distances, the particle interaction for the self-gravitating gas in physical situations is not 
gravitational. Its exact nature depends on the problem under consideration (opacity limit, Van der 
Waals forces for molecules etc.). We shall just assume a repulsive short distance potential, that is, 
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for \qi - qj\> A 
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where A << L is the short distance cut-off. 

The presence of the repulsive short-distance interaction prevents the collapse (here unphysical) of 
the self-gravitating gas. In the situations we are interested to describe (interstellar medium, galaxy 
distributions) the collapse situation is unphysical. 

2.1 The microcanonical ensemble 

The entropy of the system in can be written as 

N ■ d 3 pid 3 qi 



^,Ao=iog|^/.../n 



(2nf 



E-££--U&,...<ls) 



where 



U(qi, . . . Qn) = -Gm 2 ^ 



1 



(3) 



l<Kj<N 



M - Qj\A 



(4) 



and G is Newton's gravitational constant. 

In order to compute the integrals over the momenta pi, (1 < I < N), we introduce the variables, 



Pi = 



/2m 



Pi ■ 



We can now integrate over the angles in 3N dimensions, 
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(5) 



The delta function in the energy thus becomes the constraint of a positive kinetic energy E—Ufai, . . . qjy) > 
0. We then get for the entropy, 



\,...<t N )] 3N/2 - 1 0[E-U(&,...? N ) 



(6) 



It is convenient to introduce the dimensionless variables fj, 1 < I < N making explicit the volume 
dependence as 



qi = Lfi , n = (xi,yi,z t ) , 

< x t ,yi,zi < 1 . 

That is, in the new coordinates the gas is inside a cube of unit volume. 
The entropy then becomes 



(7) 



S(E, N) = log ■ 
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where we introduced the dimensionless variable £, 
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and 



s G m 2 N 2 
u(ri,...,r N ) = ^ 



N l<t^<N \ f l~ f 0\a 



(9) 
(10) 
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where a = A/L -C 1. 

Let us define the coordinate partition function in the microcanonical ensemble as 
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Therefore, 



S(E,N)= log 



]\f3N-2 m 9N/2-2 jJiN / '2+1 Q3N/2- 



N\r(^f) (2tt) 



37V/2 



£+-^"(ri,...,rjv) 



+ logw(£,A0 



(11) 



We can now compute the thermodynamic quantities, temperature and pressure through the standard 
thermodynamic relations 



T ~ \dE 



and p = T 



dS_ 

dV 



(12) 



where V = L stands for the volume of the system and p is the external pressure on the system. 
2.2 The canonical ensemble 

The partition function in the canonical ensemble can be written as 



N 

Z c (N,T) = ±J...Jl[ 



d 3 pi d 3 q, 



where 
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l<l<j<N 

G is Newton's gravitational constant. 

Computing the integrals over the momenta pi, (1 < I < AT) 
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We make now explicit the volume dependence introducing the variables fi, 1 < I < N defined in eq.JTJ. 
The partition function takes then the form, 



Z C (N,T) = 



1 fmTl} 



Nl V 2tt 

where we introduced the dimensionless variable r\ 

V = 

and u(fi, . . . , fV) is defined by eq.JTJJ- Recall that 

Gm 2 N 



1 N 



Y[d 3 n e V <r 1 ,...,r N ) 



1 = 1 



U 



is the potential energy of the gas. 

The free energy takes then the form, 



L 



Gm 2 N 
L T 



u(fi, . . . ,f N ) , 



(16) 



(17) 



(18) 



F = -T log Z C {N,T) = -NT log 



eV ( mT 
AT \~2^ 



3/2 



(19) 
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where 

,i ,1 n 



<Z> N (r))=log .../ UcPn , (20) 

JO JO 

The derivative of the function $^(77) will be computed by Monte Carlo simulations, mean field methods 
and, in the weak field limit 77 << 1, it will be calculated analytically. 
We get for the pressure of the gas, 

dF\ NT 77 T , 



p =-{w) T = ir-iv^). (2D 

[Here, V = L 3 stands for the volume of the box and p is the external pressure on the system.]. We see 
from ea. (|20[) that &n{v) increases with 77 since u(.) is positive. Therefore, the second term in eq.J2U is 
a negative correction to the perfect gas pressure . 

The mean value of the potential energy U can be written from eq. (|18l) as 

< U >= -T V & N (r,) (22) 

Combining eas. l|21|) and 12 I'll yields the virial theorem, 

pV_ = 1 <u> PV_ = 1 

NT 3NT ° r NT 2 3NT ' 1 ' 

where we use that the average value of the kinetic energy of the gas is ^NT. 
A more explicit form of the equation of state is 

— = 1- — &n(t)) , ( 24 ) 

NT 3N NKU ' V ' 

where 

,1 N 



&n(?i) = e ~* N[v) [ ■■■ [ jiA^i,.,4)e'" (fl - A) 
Jo Jo J=1 

1 r 1 r 1 N 1 

= -(iV-l) e - **^ / ... / T\d 3 r t - e n<?i,-,?N) _ ( 2 5) 

2 Jo Vo ^ Fi -r 2 | a 

This formula indicates that &' n (ti) is of order TV for large N. Monte Carlo simulations as well as analytic 
calculations for small rj show that this is indeed the case. In conclusion, we can write the equation of 
state of the self-gravitating gas as 

^ = m , (26) 

where the function 

m = 1 - ^ & N (v) - 

is independent of N for large N and fixed r\. [In practice, Monte Carlo simulations show that f(rj) is 
independent of N for N > 100]. 
We get in addition, 

< U >= -3NT [1 - f(rj)] . (27) 
In the dilute limit, 7/^0 and we find the perfect gas value 

/(0) = 1 • 

Equating ens. l|2"3| and (|2l)|l yields, 

$ JV (r ? )=3JV f dx 1 - . 



Relevant thermodynamic quantities can be expressed in terms of the function f(rj). We find for the free 
energy from ea. (|19f) . 

F — F ~ 3NT [ dx 1 ~ f ( X > . (28) 



10 



where 



F = -NT log 



eV (mT\ 



3/2' 



is the free energy for an ideal gas. 

We find for the total energy E, chemical potential /i and entropy S the following expressions, 



E 
S 



3NT 



dF 
ON 
dF 
df 

S Q + 3N 



T,V 



V /toT\ 3/2 1 f 



where 



is the entropy of the ideal gas. 

Notice that here the Gibbs free energy 

$ = F + pV = F + NT 



T 



/fa) -3 fdx l -J^ 



is not proportional to the chemical potential. That is, here $ ^ fiN and we have instead, 

(iN = 2NT [1-/(77)] ■ 



(29) 



(30) 



(31) 



(32) 



(33) 



This relationship differs from the customary one (see 0) due to the fact that the dilute scaling relation 
N ~ L holds here instead of the usual one N ~ L 3 . The usual relationship is only recovered in the ideal 
gas limit 77 = 0. 



The specific heat at constant volume takes the formjH], 



N \dT 



/fa) -r? /'fa) -5 



(34) 



where we used ea. <|31[) . This quantity is also related to the fluctuations of the potential energy (ALT) 2 
and it is positive defined in the canonical ensemble, 



Here, 



(AC/) 



c v = - + (Auy 



<u 2 > - <u > 2 ... , , , 

- = 3 [/fa) - 77 /' fa) - 1] 



NT 2 



The specific heat at constant pressure is given by |3| 



T 



c P = c v 



(*) 



N (dp 



and then, 



cp 



Cy + 



[/(77) " 7?/'(77)f 



/fa) + kvf'iv) 

3 4 /fa) [/fa) -77/' fa)] 

2 /fa) + |7?/'fa) 



(35) 
(36) 

(37) 



(38) 



11 



The isothermal (Kt) and adiabatic (Ks) compressibilities take the form 



K T = 





fdV 




\dp 




fdV 


V 


\~dp~ 



V 



1 



C ^K T . 
cp 



It is then convenient to introduce the compressibilities 



NT 



K 7 



1 



and ks 



NT 



K s 



cv 



v m + kvf'M V cp 

which are both of order one (intensive) in the N, L — > oo limit with N/ L fixed 
The speed of sound v s can be written as ^I! 

2 CpV 2 



kt 



(39) 



(40) 



dp 
dV 



N 



where we used ea.(l37fl in the last step. Therefore, 



T 

Nc v 



dp 
df 



dp 
dV 



T 



3 [m - Vf'(ri) - I 



(41) 



(42) 



The pressure p used in this calculation corresponds to the pressure on the surface of the system. Hence, 
this is the speed of sound on the surface of the system, this is different from the speed of sound inside 
the volume since the ground state is inhomogeneous. We compute the speed of sound as a function of 
the point in 0. 

We see that the large N limit of the self-gravitating gas is special. Energy, free energy and entropy 
are extensive magnitudes in the sense that they are proportional to the number of particles N (for 
fixed if). They all depend on the variable 77 = G ™ T N which is to be kept fixed for the thermodynamic 
limit (N — ► 00 and V — > 00) to exist. Notice that r\ contains the ratio N/L = N V^ 1 ^ which must be 
considered here an intensive variable. Here, the presence of long-range gravitational situations calls 
for this new intensive variable in the thermodynamic limit. 

In addition, all physical magnitudes can be expressed in terms of a single function of one variable: 



3 Monte Carlo Simulations 

The Metropolis algorithm^]] was applied for the first time to the self-gravitating gas in ref. 1 . We 
computed in this way the pressure, the energy, the average density, the potential energy fluctuations, the 
average particle distance and the average squared particle distance as functions of r\ for the self-gravitating 
gas in a cube of size L in the canonical ensemble at temperature T[Q . 

We implemented the Metropolis algorithm in the following way. We start from a random distribution 
of N particles in the chosen volume. We update such configuration choosing a particle at random and 
changing at random its position. We then compare the energies of the former and the new configurations. 
We use the standard Metropolis test to choose between the new and the former configurations. The 
energy of the configurations are calculated performing the exact sums as in eq. (IK)ll . We use as statistical 
weight for the Metropolis algorithm in the canonical ensemble, 

77 u(#i,...,tV) 

which appears in the coordinate partition function ea. (|20|l . The number of particles N went up to 2000. 

We introduced a small short distance cutoff A = 10~ 4 L — 10~ 8 L in the attractive Newton's potential 
according to eq. J2J) ■ All results in the gaseous phase were insensitive to the cutoff value. The partition 
function calculation turns to be much less sensible to the short distance singularities of the gravitational 
force than Newton's equations of motion for N particles. That is, solving the classical dynamics for N 
particles interacting through gravitational forces as well as solving the Boltzman equation including the 
TV-body gravitational interaction requires sophisticated algorithms to avoid excessively long computer 
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Figure 3: Average particle distribution in the gaseous phase from Monte Carlo simulations with 2000 
particles in the microcanonical ensemble for £ = —0.5, T) = 1.38, pV/[NT] — 0.277. 

times. As is clear, solving the A-body classical evolution or the kinetic equations provides the time- 
dependent dynamics and out of thermal equilibrium effects which are out of the scope of our approach. 

In the CE, two different phases show up: for r\ < tjt we have a non-perfect gas and for r\ > tjt it is 
a condensed system with negative pressure. The transition between the two phases is very sharp. This 
phase transition is associated with the Jeans instability. 

A negative pressure indicates that the free energy grows for increasing volume at constant temperature 
[see ea. <|21[) ]. Therefore, the system wants to contract sucking on the walls. 

We plot in fig. ^/(ry) = pV/[NT] as function of -q. pV/[NT] monotonically decreases with rj. 

We find an excellent agreement between the the Monte Carlo and the mean field results for pV/ [NT] 
in the gaseous phase plotted in fig. This excellent agreement appears for N as low as 500. 

In the Monte Carlo simulations the phase transition to the condensed phase happens for i] = t/t 
slightly below r\c- For N = 2000 we find rjx ~ 1.515. For r/T < i] < rjc, the gaseous phase only exists as 
a short lived metastable state as seen in our simulations. 

The average distance between particles < r > and the average squared distance between particles 
< r 2 > monotonically decrease with r\. When the gas collapses at t\t , < r > and < r 2 > exhibit a sharp 
decrease. 

The values of pVj [NT] , < r > and < r 2 > in the condensed phase are independent of the cutoff for 
a < 10~ 5 . The Monte Carlo results in this condensed phase can be approximated for -q > 2 as 



where K ~ 14. 

Since f(rf) has a jump at the transition, the Gibbs free energy $ is discontinuous and we have a phase 
transition of the zeroth order. We find from ea. (|32|l 




< r >~ 0.016 . 



(43) 



$ (collapse) - <&(?7t) 



/(collapse) - /(r? T ) ~ -21 < . 



(44) 



NT 
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We can easily compute the latent heat of the transition per particle (q) using the fact that the volume 
V stays constant. Hence, q = AE/N and we obtain from ea. (j3U|> 

q_ = Eicon***)- Efc) = g ^ ^ 2 — 3 X ?/t — —62 < . (45) 

This phase transition is different from the usual phase transitions since the two phases cannot coexist in 
equilibrium as their pressures are different. 

Ea. H43fl can be understood from the general treatment in sec. Ill as follows. We have from eas. (|24|) - 

m = 1 - 1 < - > ■ (46) 

6 r 

The Monte Carlo results indicate that < - >~ 42 is approximately constant in the collapsed region as 
well as < r > and < r 2 >. Eq. l|43|l thus follows from ea. l|4t)|) using such value of < =■ >. 

The behaviour of pV/ [NT] near rjc in the gaseous phase can be well reproduced by 

^ = f(v) "=° fc + A (47) 

where f c ^ 0.316, A ~ 0.414 and rjc ^ 1-540. 

In addition, the behaviour of (AU) 2 in the same region is well reproduced by 

(All) 2 " T =5 C C+ - ° (48) 
VVc - V 

with C ~ -1.64 and D ~ 0.901. [Notice that for finite N, (AU) 2 will be finite albeit very large at the 
phase transition]. Ea, 136|) relating f(rj) and (AU) 2 is satisfied with reasonable approximation. 

We thus find a critical region just below rjc where the energy fluctuations tend to infinity as 77 j r\c- 
The point t)t where the phase transition actually takes place in the Monte Carlo simulations is at 
r)T — 1.51 < rjc- This value for r]x is numerically close to the point where the isothermal compressibility 
kt diverges and becomes negative. As stressed in ref.^ E] these two points actually coincide. The 
collapse observed in the Monte Carlo simulations corresponds to the singularity of kt where, in addition, 
the speed of sound at the origin becomes imaginary 

Since Monte Carlo simulations are like real experiments, we conclude that the gaseous phase extends 
from 77 = till 77 = rjr in the CE and not till f] = r\c- Notice that in the literature based on the hydrostatic 
description of the self-gravitating gas |ll)l 1141 1151 , only the instability at 77 = rjc is discussed whereas 
the singularities at 77 = r\T are not considered. 



We then performed Monte Carlo calculations in the microcanonical ensemble where the coordinate 
partition function is given by ea. Hll(l . We thus used 



£+ j^u(r u ...,r N ) 



-1 3JV/2-1 



,f N ) 



as the statistical weight for the Metropolis algorithm. 

The MCE and CE Monte Carlo results coincide up to the statistical error for < 77 < 77^, that is 
for co > £ > £t — —0.19. In the MCE the gas does not clump at 77 = r\c (point C in fig. ^| and 
the specific heat becomes negative between the points C and MC. In the MCE the gas does clump 
at f ~ —0.52 , 77^ /c ~ 1.33 (point MC in fig. [JJ increasing both its temperature and pressure 
discontinuously. We find from the Monte Carlo data that the temperature increases by a factor 2.4 
whereas the pressure increases by a factor 3.6 when the gas clumps. The transition point 77^-^ in the 
Monte Carlo simulations is slightly to the right of the critical point 77MC predicted by mean field theory. 
The mean field yields for the sphere 77MC = 1-2598 . . .. 



As is clear, the domain between C and MC cannot be reached in the CE since cy > in the CE as 
shown by eg. 1 (55)1 . 

We find an excellent agreement between the Monte Carlo and Mean Field (MF) results (both in the 
MCE and CE). (This happens although the geometry for the MC calculation is cubic while it is spherical 
for the MF). The points where the collapse phase transition occurs (rjr and i]J IC ) slowly increase with 
the number of particles N. 
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Canonical, Gas Phase, eta = 1 .5 
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Figure 4: Average particle distribution in the gaseous phase from Monte Carlo simulations in the canonical 
ensemble for T) = 1.5 and ./V = 2000 



We verified that the Monte Carlo results in the gaseous phase (r) < rjx) are cutoff independent for 
10" 3 > a > IQ- 7 . 

As for the CE, the Gibbs free energy is discontinuous at the transition in the MCE. The transition is 
then of the zeroth order. We find from eg. 1321) 

^collapse) = /( ^ ^ Q ? > Q 

-< * ^ gas -L gas 

where we used the numerical values from the Monte Carlo simulations. Notice that the Gibbs free energy 
increases at the MC transition whereas it decreases at the C transition [see ea. (!44ll ]. 

Here again the two phases cannot coexist in equilibrium since their pressures and temperatures are 
different. 

We display in figs. 13121 the average particle distribution from Monte Carlo simulations with 2000 
particles in the microcanonical ensemble at both sides of the gravothermal catastrophe, i. e. 77 = t\mc- 
Fig. |3 corresponds to the gaseous phase and fig. |3to the collapsed phase. The inhomogeneous particle 
distribution is clear in fig. [3] whereas fig. |21 shows a dense collapsed core surrounded by a halo of particles. 

The different nature of the collapse in the CE and in the MCE can be explained using the virial 
theorem [see ea. (|23f) ] 

pv =l , U 



NT NT 

When the gas collapses in the CE the particles get very close and U becomes large and negative while T 
is fixed. Therefore, jpp may become large and negative as it does. 
We can write the virial theorem also as, 

p v-\nt = \e. 
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Canonical, Collapsed Phase eta = 1 .53 



+ one particle 




Figure 5: Average particle distribution in the collapsed phase from Monte Carlo simulations with 2000 
particles in the canonical ensemble for -q — 1.53, pV/[NT] = —14.44. 

When the gas is near the point MC, E < is fixed and we have T > 0. Therefore, jpp as well as 
U = E — 3 NT/2 cannot become large and negative as in the CE collapse. This prevents the distance 
between the particles to decrease. Actually, the Monte Carlo simulations show that < r > increases by 
18% when the gas collapses in the MCE. 

Figs. 0]andO3depict the average particle distribution from Monte Carlo simulations with 2000 particles 
in the canonical ensemble at both sides of the collapse critical point, i. e. rj = r\c- Fig. 01 corresponds to 
the gaseous phase and fig. EJto the collapsed phase. The inhomogeneous particle distribution is clear in 
fig. 01 whereas fig. [S] shows a dense collapsed core surrounded by a very little halo of particles. 

Notice that the collapsed phases are of different nature in the CE and MCE. The core is much tighter 
and the halo much smaller in the CE than in the MCE. 

Figs. 3 and 5 depict the average particle distribution for the gaseous phase in the MCE and the CE, 
respectively. In this phase, the MC simulations give identical descriptions for large N in both ensembles. 
[This important point will be further demostrated in sec. VI by functional integral methods]. The average 
configurations in figs. 3 and 5 describe a self-gravitating gas in thermal equilibrium within a cube. We 
may call it the isothermal cube by analogy with the well known isothermal sphcre[TT ] -|16 | . 

4 Mean Field Approach 

Both in the microcanonical and the canonical ensembles the coordinate partition functions are given by 
3A-uple integrals [eas. (|llfl and l|20|) . respectively]. In the N — > oo limit both 3A-uple integrals can be 
recasted as functional integrals over the continuous particle density as shown in refs.[fl|2|. We present 
here the mean field approach in the canonical ensemble. The mean field approach for the microcanonical 
and grand canonical ensembles can be find in refs.[D|2]. 

The coordinate partition function e* N ('') given by ea. (|20|) can be recasted as a functional integral in 
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the thermodynamic limit. 

e **(„) N^>i I I Dpd&e -N.aW.),i,v] (49) 



fj f x d^y f ( f 

sc\p{.),a,rj\ = -- J |^_^-| P(x)P(y)+ J d 3 xp{x) log p{x) - id f J d x p(x) - 1 

where we used the coordinates a? in the unit volume. The first term is the potential energy, the second 
term is the functional integration measure for this case (see appendix A). Here N p(x) stands for the 
density of particles. 

The integration over a enforces the number of particles to be exactly N: 

d 3 x p(x) = 1 (50) 

That is, in the coordinates q (running from to L), the density of particles is 



N 
L3 



r TV 

p(q) with / d 3 q — p(q) = N . 



The functional integral in ea. (|49|l is dominated for large N by the extrema of the 'effective action' 
sc[p(-), &, fj], that is, the solutions of the stationary point equation 

i r->\ f d3 V Ps(y) f xt\ 
\ogp s {x) - J) \ — — — = a s , (51) 

J \x-y\ 

a = id is a Lagrange multiplier enforcing the constraint (|5UI) . 
Applying the Laplacian and setting <f>(x) = \ogp s (x) yields, 

X7 2 <f>(x)+4nr]e' l > (S) =0, (52) 

This equation is scale covariant j^j. That is, if 4>(x) is a solution of ea- H52|l . then 

A (f) = 0(Ax)+logA 2 (53) 

where A is an arbitrary constant is also a solution of ea. (|52J) . For spherically symmetric solutions this 
property can be found in ref.[H]. 

Integrating eq.(jS2J over the unit volume and using the constraint 1)50(1 yields 

V(j)(x) ■ ds = (54) 

where the surface integral is over the boundary of the unit volume. 

In the mean field approximation we only keep the dominant order for large N. Therefore, only the 
exponent at the saddle point accounts and according to ea. H19|) we find for the free energy 

F = F + NTs(ri) + O(N ) 

Hence, in the mean field approximation, the function f(rj) is given by 

/m^)-1 + ||, (56) 

From ea. (|49fl we can compute s(rj) in terms of the saddle point solution as follows 

Tj f d^x d^y f 
s(v) = sc[p s (-),a s ,rj\ = -~ J p s (x) p s {y) + J d 3 x p s {x) \ogp s {x) . (57) 

Using ea. l|51|) we find the equivalent expression, 

s(v) = y + \ f m d 3 x . (58) 
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Figure 6: fMF(i] R ) — PV/[NT] as a function of rj R in the MF approximation [ea. (|64[l ]. fMF(i] R ) 
has a square root branch point at i] R . The points GC, C and MC indicate the transition to the col- 
lapsed phase for each ensemble (grand canonical, canonical and microcanonical, respectively): rj R c = 
0.797375 . . . , r] R = 2.517551 . . . , rj^c — 2.03085 . . . (notice that ij^j C is in the second Riemann sheet). 
Since E/[3NT] = fMF(v R ) ~ \ j this plot also shows the energy per particle as a function of r\ R . Further- 
more, the particle density at the surface is proportional to /m_f('7 R ) [see ea. H62[) ]. 



5 Mean Field Results 

Let us summarize here the main results of ref. in the mean field approach for spherical symmetry. 
The saddle point is given by 



(r) = logp(r) = log 



A 



Here p{r) is the particle density and x(A) obeys the equation 



+ X (Ar) 



x"(A) + x x'(A)+e*( A )=0 , X '(0)=0 

x{x) is independent of rj , and A is related to rj R through 

A x'(A) = -V R ■ 

We have in addition, 

0(1) = log 



X(0) = ; . 



3/mf(^) 



47T 



p(l) = — f MF (v R ) 

47T 



(59) 

(60) 

(61) 
(62) 



where 



fMF(V R ) 



3n R 



The function fMF(r] R ) obeys the Abel equation, 



(63) 



ti r (2>!mf - l)f M F(v R ) + (3/mf - 3 + V R )f MF = 



(64) 



18 




19 



For r] R = it follows from eg. 1(64(1 that 

/mf(0) = 1 . (65) 

Integrating ea. 1(64(1 with respect to r/ R yields, 

3 / — [1 - !mf{x)\ = 3[f MF (v R ) - 1] + V R - logf M F(v R ) 
Jo x 

We derive in £Q the properties of the function fMF(v R ) from the differential equation 1(64(1. One easily 
obtains for small rj R (dilute regime), 

fM F (v R ) = i- 1 f~ i -^+o([ V R r). 

These terms exactly coincide with the perturbative calculation in the dilute regime for spherical symmetry 

m 

We plot in fig. ^Jmf(v R ) as a function of rj R obtained by solving ea. 1(64(1 by the Runge-Kutta method. 
We see that /mf(^?' R ) is a monotonically decreasing function of i] R for < r\ < i] R . At the point 
i] R = rj R , the derivative f'Mp(T) R ) takes the value — oo. It then follows from ea. ((64(l that 

fMFiVc) = \ ■ 

At the point rj R the series expansion for fMF(v R ) m powers of r/ R diverges. Both, from the ratio test on 
its coefficients and from the Runge-Kutta solution, we find that 

77^ = 2.517551... (66) 

From ea. (|64|l we find that fMF(v R ) ~ | has a square root behaviour around rj R = r] R : 



Inserting the numerical value 1)66(1 for 77^ yields, 



f M F(v R ) ' ,R = S I + 0.213738 ...^J V R - V R + 0.172225 . . . fag - r, fi ) + O |fa« - ,7«) 3 / 2 j (67) 

We see that fMF{n R ) becomes complex for r] R > rj R . Recall that in the Monte Carlo simulations the gas 
phase collapses at the point rj R < rj R . fMF{n R ) is a multivalued function of ij R as well as all physical 
magnitudes [see eq.JBSJ]. 

The points GC, C and MC correspond to the collapse phase transition in the grand canonical, canonical 
and microcanonical ensembles, respectively. Their positions are determined by the breakdown of the mean 
field approximation through the analysis of the small fluctuations. 

As noticed before, the CE only describes the region between the ideal gas point, r] R = and C in fig. 
1. The MCE goes beyond the point C (till the point MC) with the physical magnitudes described by 
the second sheet of the square root in eas. (|67|l (minus sign). We have near C between C and MC, 

f MF (v R ) "*= S I - 0-213738 . . . y/ng - V R + 0-172225 . . . ( V R - V R ) + O [(rj R - V R f/ 2 

The function fMF(v R ) takes its absolute minimum at r] R = ij R in = 2.20731 ... in the second sheet 
where fMF(i] R in) — 0.264230.... Since fMF{f] R ) < I implies that the total energy is negative [see 
ca. (168fi ] . the gas is in a 'bounded state' for r] R beyond r\ R = 2.18348 ... in the first sheet. 



For the main physical magnitudes in the mean field approach, [that is, from the above saddle point 
and neglecting the fluctuations around it] we find: 

P v f ( Rs 

Jjj, = JMF(V ) 

= Z[l-fMF(V R )]-V R + l0gf M F(V R ) 
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Figure 9: (cv)mf as a function of r] R from mean field ea. lj7U|) . Notice that (cv)mf diverges at the point 
C, that is for r/ R = 2.517551 . . . 



S-Sg 
N 



= 6[f M F{v R ) ~ 1] +ri R - log f M F(v R ) 



(68) 



= 3[/mf (if ) - \] ■ 
We find for the speed of sound squared at the surface 



T 



fMF(v R ) 



3/mf(^) + % 



11 s 

2 / 



f M F(v R ) + h 



The specific heat at constant volume takes the form 



(cv)mf = 6 fuF (v R ) ~l + V R + ~ f —, £ ; 

2 3 Jmf(V ) - 1 



(69) 



(70) 



We plot in Fig. I§1 ea. (|70|l for (cv)mf as a function of r\. We see that (cv)mf increases with r] till it 
tends to +oo for r] R | ?7c- ^ nas a square-root branch point at the point C. In the stretch C-MC (only 
physically realized in the microcanonical ensemble), (cv)mf becomes negative. We shall not discuss here 
the peculiar properties of systems with negative Cy as they can be find in refs. |12l IT!T1 116| 

From eas. (|67[) and (|Tf)|> we obtain the following behaviour near the point C in the positive (first) 
branch 



(c v )mf n = ° 0.80714. 



, (jf - rfy 1 ! 2 - 0.19924 . . . + 0(^r]§ - r) R ) 



(71) 



and between C and MC in the negative (second) branch 

(cv)mf -0.80714 . . . fag - v R y 1/2 ~ 0.19924 . . . + 0(^ V R - t?«) 

Finally, (cv)mf vanishes at the point MC r] R lc = 2.03085 . . .. 

The isothermal compressibility in mean field follows from eas. (|39ll and (|64() 



(«t) 



MF 



2f MF (v R ) 



r, R - 2 



6f M F(v R ) -V 1 



(72) 
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Figure 10: (kt)mf as a function of 77 from mean field ea. (|72[l . Notice that (kt)mf diverges at r] R 
r] R = 2.43450 . . . 



We plot (kt)mf in fig. HOI We see that (k t )mf is positive for < rj R < 77^ = 2.43450 . 
diverges. The point rj R is defined by the equation 



6fMF(Vr) 



. 



where (kt)mf 



(73) 



We find from eqs.lHO} and (J73J) that 
{kt)mf diverges for rj ~ rj R as 



/mf(1t) — o 



(74) 



(«t) 



mf 



n^n* 9 (77I - 2) 



4n§(ri$ -r) R ) 



•0(1) 



0.40157. . 
r/y — 77^ 



0(1) 



(kt)mf is negative for 77^ < r7 fl < 77^ and exactly vanishes at the point C. (kt)mf then becomes 
positive in the stretch between C and MC only physically realized in the microcanonical ensemble. 

Notice that the singularity of {kt)mf at 77^ = 77I? = 2.43450. . . is before but near the point C. It 
appears as a preliminary signal of the phase transition at C . rj R is the transition point 77T seen with the 
Monte Carlo simulations (see fig. 1). (Recall that 77^ ~ 1.515 corresponds to 77^ ~ 2.44). 



It is easy to understand the meaning of a large compressibility. From the definition (|39|) 

SV . VSp 

— = -K T 6p = - KT — . 



(75) 



A large compressibility implies that a small increase in the pressure (dp <C NT/V) produces a large 
change in the density of the gas. That means a very soft fluid. 

For negative compressibility, eq.J75J) tells us that the gas increases its volume when the external 
pressure on it increases. This is clearly an unusual behaviour that leads to instabilities as one sees by 
computing the velocity of sound as a function of r^E]. 

For the specific heat at constant pressure we find 



(cp)mf = l2f MF ( V R ) 



3 24 {v R - 2) f MF (v R ) 
2 



6 Imp (v r ) 



77 



R 



(76) 
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Defining Equation 




PHYSICAL MEANING 


GC 


1.7772... 


0.797375 . . . 


2-3r/g c f MF ( v K c )=0 


0.836076... 


Collapse in the GCE. 


3 


3.38626 . . . 


1.73745... 


3- V R + x(A) = 


0.622424... 


Energy density 
vanishes at r = 0. 


2 


4.73739 . . . 


2.18348... 


2 !mf{^)- 1 = 


0.5 


Total Energy vanishes. 


T 


6.45077... 


2.43450 . . . 


G fMF(rtf) - rtf = 


0.40575 . . . 


kt and cp diverge. 
Collapse in the CE. 


C 


8.993195... 


2.517551... 


1 - 3 f MF (Vc) = 


1/3 


cy diverges. 


Min 


22.5442... 


2.20731... 


f M F^WmiTi) = 


0.264230... 


Minimum of 

P V/[NT] 

in the gas phase 


1 


25.7991... 


2.14675... 


^P MF (Vi ) - (38 - 8/7f ) jW(??f )+ 
+r?f = 


0.265290... 


Vg and cp vanish. 


MC 


34.36361 . . . 


2.03085 . . . 


12f 2 MF (viic) - (11 - 2< c ) /mf(4)+ 
+1 = 


0.273512... 


Collapse in the MCE. 
cy vanishes. 
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Figure 11: Mean value of the distance between particles (< r >) and mean value of the squared distance 
between particles (< r 2 >) as functions of rj R in the mean field approach from eas. H79(l . Notice that the 
particles are inside a sphere of unit radius. 



These average distances are defined as 



< r > 



< r" > = 



'-f'\ <p{r)p(r')> d 3 r d 3 



'- < p(r) p(f') > d 3 r d 3 r' 



In the mean field approximation we have 



< p(r) P(r") >= Pmf{t) Pmf{t') + {j-j 



(77) 



(78) 



In addition, in the spherically symmetric case we use ea. (|80|l for the particle density. In Appendix E we 
compute the integrals in eas. l77f) and we get as result, 



<r > = 2-2 



r dr 



(r)-^(l) 



< r 2 > 



12 



2 _ — / r* dr [0(r) - 0(1)] 
V Jo 



(79) 



where <f>(r) is given by ea. l|59|) . 

We plot < r > and < r 2 > as functions of ri R in fig. ^] Both < r > and < r 2 > monotonically 
decrease with X(r] R ). Their values for the ideal gas are 



< r >\ v=0 - | = 1.02857 



< r > 



r)=0 



At the critical points (C for the canonical ensemble and MC for the microcanonical ensemble) the average 
distances sharply decrease. Both < r > and < r 2 > have infinite slope as functions of r/ R at the point C. 

We plot in fig. ^Jthe Monte Carlo results for < r > in a unit cube together with the MF results in a 
unit sphere. Notice that < r > sharply falls at the point T clearly indicating the transition to collapse. 
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Figure 12: Mean value of the distance between particles (< r >) in a unit cube from Monte Carlo 
simulations and in a unit sphere from mean field as functions of rj R . 



5.2 Particle Distribution 

The particle distribution at thermal equilibrium obtained through the Monte Carlo simulations and mean 
field methods is inhomogeneous both in the gaseous and condensed phases. In the dilute regime i)< 1 
the gas density is uniform, as expected. 

For the Monte Carlo simulations the density of particles in the cube are shown for the gaseous and 
for the condensed phases in figs. 3-6, respectively. 

In the mean field approximation and for the spherically symmetric case, the particle density is given 

by 

j., ^ A 2 e*( Ar ) 

Pmf(t) = e+W = - A =- , 0<r<l. (80) 



The mass inside a radius r is then given by 



M(r) = 4ir f r' 2 p MF (r') dr' = -^- x '(Ar) 
Jo V 



where we used ea. (|6()|l . For small r this gives 



(81) 



We find an uniform mass distribution near the origin. This is simply explained by the absence of gravi- 
tational forces at r = 0. Due to the spherically symmetry, the gravitational field exactly vanishes at the 
origin. The particles exhibit a perfect gas distribution in the vicinity of r = 0. Actually, ea. (|81|l is both 
a short distance and a weak coupling expression. Ea. (|81|l is valid in the dilute limit i] R <C 1 for all 
< r < 1. 
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Figure 13: log M(r) vs. log r for five different values of rj R : i] R IC — 2.03085 . . . , rj R = 2.517551 . . . , rj R = 
2, r?g c = 0.797375 ... and t] R = 0.1. 





D 


C 


0.1 


2.97 


1.0 


n% c = 0.797375 . . . 


2.75 


1.03 


2.0 


2.22 


1.1 


7$ = 2.517551 . . . 


1.60 


1.07 


rj R c = 2.03085 . . . 


0.98 


1.11 



TABLE 2. The Fractal Dimension D and the proportionality coefficient C as a function of 7y fi from 
a fit to the mean field results according to M. (r) ~ C r D . 

We plot in fig. ^|the particle distribution for 90% of the particles for several values of rj R . We exclude 
in the plots the region M{r) < 0.1 where the distribution is uniform. 

We find that these mass distributions approximately follow the power law 

M(r)~Cr D (82) 

where, as depicted in Table 2, D slowly decreases with X(rj R ) from the value D = 3 for the ideal gas 
(i] = 0) till D = 0.98 in the extreme limit of the MC point. 



5.3 The speed of sound as a function of r 

For very short wavelengths A s <C L, the sound waves just feel the local equation of state Q and the 
speed of sound will be that of an ideal gas. For long wavelengths (of the order L), the situation changes. 
The calculation in eq. (|69|l corresponds to the speed of sound for an external wave arriving on the sphere 
in the long wavelength limit. Let us now make the analogous calculation for a wave reaching the point q 
inside the gas. 

Our starting point is eq.JHJ, 

2 c P V 2 fdp(q)\ 

vM = -^N\-mr) T ^ (83) 
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where cp and cy are the specific heats of the whole system at constant (external) pressure and volume, 
respectively, and p(q) is the local pressure at the point q. The local pressure in the spherically symmetric 
case can be written in a more explicit way using eqs.JSSJ) and 



p(r) V 



NT 3t]' 

We find for the spherically symmetrical case in MF 
i> 2 (r) c P A 2 



(84) 



T c v 9rj R [3f(r] R ) - 1] 
where we used eas.l(8ijl. and 



6 f(r) R ) + Xr x'(A r)] e^ x r > , (85) 



dr]\ r\ 
dV) T ~ ' 

[A is a function of i] R as defined by ea. ffilfl ]. 

At the surface, (r = 1), v^(r) reduces to eq. (f^|) after using eqs.lJBTJj (j7U|) and f75|l. 

For r] R < r\ R = 2.43450 ... in the first sheet, v 2 (r) is positive and decreases with r as shown in fig. 

At r\ R = rj R , v 2 s (r) diverges for all < r < 1 due to the factor cp in eq.JSSJ [cfr. ea. l|76|l ]. The 
derivative of p with respect to V is proportional to 6 f(rj R ) + Xr x'(Ar) as we see in ea. l|85|) . At r = 1 
this factor becomes 6 f{r] R ) — rj R [see ea. l|61|l ] which exactly vanishes at r/ R — r/ R canceling the singularity 
that cp possesses at such point [see ea. l(7T>|) ]. Thus, w 2 (l) is regular at r] R = r] R . 

v1(r) becomes large and positive below and near r] R and large and negative above and near 
r] R = 2.43450 ... as depicted in fig. This singular behaviour witness the appearance of strong 

instabilities at r/ R — r/ R . For r] R > r\ R 1 v s (r) becomes imaginary indicating the exponential growth of 
disturbances in the gas. This phenomenon is especially dramatic in the denser regions (the core). 

For r/ R beyond r/ R and before r] R = 2.14675 ... in the second sheet w 2 (r) stays negative around the 
core while it becomes positive in the external regions For example, f 2 (r) is positive at r/ R — rj R for 
r > 0.4526.... 

For rj R in the second sheet beyond r] R and before rj R iCl f 2 (r) is positive in the core and negative 
outside. 
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Figure 15: The speed of sound v 2 (r)/v 2 (l) as a function of r for values of r/ R around r/ R — 2.43450 . . .. 
Positive values of v 2 (r) / v 2 (1) correspond to rj R < r\ R and negative values of v 2 (r)/v 2 (l) correspond to 
T\ R > rj R . We see that the speed of sound squared tends to +oo in the bulk (r < 1) for r] R | rj R while it 
tends to — oo for rj R { rj R . 



6 ^-dimensional generalization 

The self-gravitating gas can be studied in ^-dimensional space where the Hamiltonian takes the form 



N 



1=1 



El. 

2m 



Gm 2 



E 



and 



^-Gm 2 V log— — , ati/ = 2. 

2m \qi — qAa 

1=1 l<Kj<N |y( HnA 



(86) 



(87) 



The partition function in the microcanonical, canonical and grand canonical ensembles takes forms 
analogous to eqs.(II.2), (III. 1) and (VI. 7) in paper I, respectively. 
We now find for the microcanonical ensemble, 



S(E,N)= log 



jyvN-2 m 3vN/2-2 jv(2-v /2) N+v-2 QvN/2-1 



where the coordinate partition function takes now the form 



+ logw{£,N) . 



r l r l N 



— u(fi,...,rjv) 



vN/2-1 



t. + jju{r 1 ,...,r N ) 



with 



and 



EL 



v-2 



Gm 2 N 2 
u{r u ...,r N ) = ^- 



(88) 



l<l<j<N 



\n-Ti\a 



v-2 ■ 



28 




-0.8 1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 

Figure 16: The S-wave determinants Dgc(v R )i ^c(v R ) an( i Dmc(v R ) m the grand canonical, canonical 
and microcanonical ensembles, respectively, as functions of r] R . Notice that the mean field approximation 
for each ensemble breaks down as soon as the respective determinant becomes negative when r\ increases 
starting from 77 = 0. 



In the canonical ensemble we obtain now, 



Nl V 2tt 
where the variable 77 takes the form 







r[ e n «(fi,...,fjv) 



1=1 



V = • (89) 



Gm z N 
L v - 2 T 

As we can see from eas. (|88|l - (189t the only change going off three dimensional space is in the exponent of 
L. 



In ^-dimensional space the thermodynamic limit is defined as V, N — > 00 keeping r\ and £ fixed. 
That is, N/L v ~ 2 = N/V 1 " 2 /" is kept fixed. The volume density of particles N/V vanishes as V~ 2 I V for 
V, N — > 00 and v > 2. It is a dilute limit for v > 2. 

When v < 2, one has to assume that the temperature tends to infinity in the thermodynamic limit in 
order to keep 77 and £ fixed as V, N — > 00. 

7 Functional Determinants: the Validity of Mean Field 

The mean field gives the dominant behaviour for N — > 00. The Gaussian functional integral of small 
fluctuations around the stationary points was evaluated in ref.pj]. 

As remarked in pQ, the three statistical ensembles (grand canonical, canonical and microcanonical) 
yield identical results for the saddle point. However, the small fluctuations around the saddle take 
different forms in each ensemble. 

Adding the contributions from the functional determinant to the mean field results ea. (|68J) in the 
grand canonical, canonical and microcanonical ensembles corresponds to include 1 /N corrections in the 
function f(r) R ) as follows, 

n R d / 1 \ 

f(il R ) = fMF(v R ) + \°g D Gc(v R ) + O T?2 grand canonical ensemble; 



6Nd V R b lJ ° w ' \N 2 
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/(r/ K ) = /MHf) + f^^R ^SDC(V H ) + O yjpj canonical ensemble; (90) 

ri R d ( 1 \ 

/07 ) = !mf(j] R ) + lo S d mc{v ) + C ( jy2 ) microcanonical ensemble. 

where Dec, Dc and Dmc are the functional determinants in the grand canonical, canonical and micro- 
canonical ensembles, respectively. We plot in fig. ^J|the S-wave part of these determinants as functions 
of rj R . The higher partial waves give a positive definite contribution to these three determinants. 
We find in the grand canonical ensemble that 

r, R -^ logDetcofa*) -^c_ ^ 

Therefore, the energy and the entropy tend to minus infinity when r) R f r/g C . This behaviour 
correctly suggests that the gas collapses for rj R f Indeed, the Monte Carlo simulations yield a large 

and negative value for in the collapsed phase [I]. 

We want to stress that the mean field values provide excellent approximations as long as N\r] R c — 
rj R \ >> 1 in the grand canonical ensemble. Namely, the mean field is completely reliable for large N 
unless r) R gets at a distance of the order A" 1 from r] R c . 

The clumping phase transition in the canonical ensemble takes place when Dc{ij R ) vanishes at r) R = 
r) R . Near such point the expansion in 1/A breaks down since the correction terms in ea. (|90J) become 
large. Mean field applies when N\r) R — rj R \ » 1. 

Since 

\ *F** Dc ^ 1 ^ "°° • (91) 

cci. (I91|I correctly suggests that PV/[NT], E/[NT) and the entropy per particle become large and 
negative for rj R | rj R . Indeed the Monte Carlo simulations yield a large and negative value for these 
three quantities in the collapsed phase 1 . 
For i] R i 

Vmc reaching the point MC, we find 



6 dr, R b LW 7 12{ V R - V R C ) 

We see that the MF predicts that pV/[NT] grows approaching the critical point MC. This behaviour is 
confirmed by the Monte Carlo simulations. At the point MC, pV/[NT] increases discontinuously by 50% 
in the Monte Carlo simulations. 



8 The Interstellar Medium 

The interstellar medium (ISM) is a gas essentially formed by atomic (HI) and molecular (H2) hydrogen, 
distributed in cold (T ~ 5 — 50-ftT) clouds, in a very inhomogeneous and fragmented structure. These 
clouds are confined in the galactic plane and in particular along the spiral arms. They are distributed in a 
hierarchy of structures, of observed masses from 10 -2 Mq to 10 6 M Q . The morphology and kinematics of 
these structures are traced by radio astronomical observations of the HI hyper fine line at the wavelength of 
21cm, and of the rotational lines of the CO molecule (the fundamental line being at 2.6mm in wavelength), 
and many other less abundant molecules. Structures have been measured directly in emission from O.Olpc 
to lOOpc, and there is some evidence in VLBI (very long based interferometry) HI absorption of structures 
as low as 10~ 4 pc — 20 AU (3 10 14 cm). The mean density of structures is roughly inversely proportional 
to their sizes, and vary between 10 and 10 5 atoms /cm 3 (significantly above the mean density of the 
ISM which is about 0.1 atoms /cm 3 or 1.6 10~ 25 g/cm 3 ). Observations of the ISM revealed remarkable 
relations between the mass, the radius and velocity dispersion of the various regions, as first noticed by 
Larson (221 ^ an d since then confirmed by many other independent observations (see for example ref.pS)). 
From a compilation of well established samples of data for many different types of molecular clouds of 
maximum linear dimension (size) R, total mass M and internal velocity dispersion Av in each region: 

M(R) - R dH , Av ~ R q , (92) 

over a large range of cloud sizes, with 10~ 4 — 10~ 2 pc < R < 100 pc, 

1.4 < d H < 2, 0.3 < q < 0.6 . (93) 
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These scaling relations indicate a hierarchical structure for the molecular clouds which is independent 
of the scale over the above cited range; above 100 pc in size, corresponding to giant molecular clouds, 
larger structures will be destroyed by galactic shear. 

These relations appear to be universal, the exponents dn , q are almost constant over all scales of the 
Galaxy, and whatever be the observed molecule or element. These properties of interstellar cold gas are 
supported first at all from observations (and for many different tracers of cloud structures: dark globules 
using 13 CO, since the more abundant isotopic species 12 CO is highly optically thick, dark cloud cores 
using HCN or CS as density tracers, giant molecular clouds using 12 CO, HI to trace more diffuse gas, 
and even cold dust emission in the far- infrared) . Nearby molecular clouds are observed to be fragmented 
and self-similar in projection over a range of scales and densities of at least 10 4 , and perhaps up to 10 6 . 

The physical origin as well as the interpretation of the scaling relations l|92|l have been the subject of 
many proposals. It is not our aim here to account for all the proposed models of the ISM and we refer 
the reader to refs.["""3] for a review. 

The physics of the ISM is complex, especially when we consider the violent perturbations brought by 
star formation. Energy is then poured into the ISM either mechanically through supernovae explosions, 
stellar winds, bipolar gas flows, etc.. or radiatively through star light, heating or ionizing the medium, 
directly or through heated dust. Relative velocities between the various fragments of the ISM exceed 
their internal thermal speeds, shock fronts develop and are highly dissipative; radiative cooling is very 
efficient, so that globally the ISM might be considered isothermal on large-scales. Whatever the diversity 
of the processes, the universality of the scaling relations suggests a common mechanism underlying the 
physics. 

We proposed that self-gravity is the main force at the origin of the structures, that can be perturbed 
locally by heating sources 01 EJ. Observations are compatible with virialised structures at all scales. 
Moreover, it has been suggested that the molecular clouds ensemble is in isothermal equilibrium with 
the cosmic background radiation at T ~ 3K in the outer parts of galaxies, devoid of any star and 
heating sources. This colder isothermal medium might represent the ideal frame to understand the role 
of self-gravity in shaping the hierarchical structures. 

In order to compare the properties of the self-gravitating gas with the ISM it is convenient to express 
m, T and L in 77 in appropriate units. We find from ea. i|17|) 

„ = 0.52193 ^ , 

where m is in multiples of the hydrogen atom mass, T in Kelvin, L in parsecs and Mq is the mass of 
the cloud in units of solar masses. Notice that L is many times (~ 10) the size of the cloud. 

The observed parameters of the ISM clouds [2*5] yield an 77 around ~ 2.0 for clouds not too large: 
M.® < 1000. Such 1) is in the range where the self-gravitating gas exhibits scaling behaviour. 

We conclude that the self-gravitating gas in thermal equilibrium well describe the observed fractal 
structures and the scaling relations in the ISM clouds [see, for example fig. 1131 and table 2]. Hence, 
self-gravity accounts for the structures in the ISM. 



9 Discussion and Conclusions 

We have presented here a set of new results for the self-gravitating thermal gas obtained by Monte Carlo 
and analytic methods. They provide a complete picture for the thermal self-gravitating gas. 

Contrary to the usual hydrostatic treatments we do not assume here an equation of state 

but we obtain the equation of state from the partition function [see eq. (J2J5J)] . We find at the same time 
that the relevant variable is here rj R = Gm 2 N /\V 1 ^T]. The relevance of the ratio Gm 2 /[1/ 1 / 3 T] has 
been noticed on dimensionality grounds |10j . However, dimensionality arguments alone cannot single out 
the crucial factor N in the variable rj R . 

The crucial point is that the thermodynamic limit exist if we let N — > 00 and V — > 00 keeping r) R 
fixed. Notice that 77 contains the ratio N I/ -1 / 3 and not N/V. This means that in this thermodynamic 
limit V grows as N 3 and thus the volume density p = N/V decreases as ~ N~ 2 . 77 is to be kept fixed for a 
thermodynamic limit to exist in the same way as the temperature. pV, the energy E, the free energy, the 
entropy are functions of 77 and T times N. The chemical potential, specific heat, etc. are just functions 
of 77 and T. 

Starting from the partition function of the self-gravitating gas, we have proved from a microscopic 
calculation that the local equation of state p(r) — T pv{r) an d the hydrostatic description are exact. 
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Indeed, the dilute nature of the thermodynamic limit (N ~ L — ► oo with N/L fixed) together with the 
long range nature of the gravitational forces play a crucial role in the obtention of such ideal gas equation 
of state. 

We find collapse phase transitions both in the canonical and in the microcanonical ensembles. They 
take place at different values of the thermodynamic variables and are of different nature. In the CE the 
pressure becomes large and negative in the collapsed phase. The phase transition in the MCE is sometimes 
called 'gravothermal catastrophe'. We find that the temperature and pressure increase discontinuously 
at the MCE transition. Both are zeroth order phase transitions (the Gibbs free energy is discontinuous). 
The two phases cannot coexist in equilibrium since the pressure has different values at each phase. 

The parameter rj R [introduced in ea. i|17Jl ] can be related to the Jeans length of the system 

fW 1 , , 

dj = \ /==, 94 

V m \JG m p 

where p = N/V stands for the number volume density. Combining eas. l|17(l and l|94[l yields 




We see that the phase transition in the canonical ensemble takes place for dj ~ L. [The precise numerical 
value of the proportionality coefficient depends on the geometry] . For d j > L we find the gaseous phase 
and for dj < L the system condenses as expected. Hence, the collapse phase transition in the canonical 
ensemble is related to the Jeans instability. 

The latent heat of the transition (q) is negative in the CE transition indicating that the gas releases 
heat when it collapses [see ea. (|45H ]. The MCE transition exhibits an opposite behaviour. The Gibbs free 
energy increases at the MCE collapse phase transition (point MC in figQ) whereas it decreases at the CE 
transition [point T in fig. see ea. (|44|l ]. Also, the average distance between particles increases at the 
MCE phase transition whereas it decreases dramatically in the CE phase transition. These differences 
are related to the MCE constraint keeping the energy fixed whereas in the CE the system exchanges 
energy with an external heat bath keeping fixed its temperature. The constant energy constraint in the 
MCE keeps the gas stable in a wider domain and makes the collapse transition softer than in the CE. 
Notice that the core is much tighter and the halo much smaller in the CE than in the MCE [see figs. |21 
and E]. 



9.1 Hydrodynamical description 

More generally, one can investigate whether a hydrodynamical description will apply for a self-gravitating 
gas. One has then to estimate the mean free path (I) for the particles and compare it with the relevant 
scales a in the system [25]. We have, 

i~ — (95) 
p v <Jt jy &t 

where py = -r^ p is the volume density of particles and at the total transport cross section. 

Due to the long range nature of the gravitational force, a t diverges logarithmically for small angles. 
On a finite volume the impact parameter is bounded by L and the smaller scattering angle is of the order 
of 

Aq G m 2 
~q LT~ 

since q = mv ~ VmT and Aq ~ - TvP' 

We then have for the transport cross section[2*5). 

(Gm) 2 , LT ( L N\ 2 N 
\v — 1/| 4 G m z \ T) J 77 

where we used that v ~ an d \v — tf\ ~ Aq/m. As we see, the collisions with very large impact 

parameters (~ L) dominate the cross-section. 

From eas. ()95|l and (|96f) . we find for the mean free path: 



L ( Gm 2 \~ 1 L rj 



2 



l ~ N\TL) i0 g^~ AT3log^ ' (97) 
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where we have here replaced pv by -p-. A more accurate estimate introduces the factor p = in the 
denominator of I. This factor for spherical symmetry can vary up to two orders of magnitude [see 
but it does not change essentially the estimate lj9*7)l. 

We see from ea. (|97fl that in the thermodynamic limit I becomes extremely small compared with any 
length a = O(N ) that stays fixed for N — > oo. We find from eq.JHU), 

I 1 

a~ A^loif ' 

In conclusion, the smallness of the ratio I /a (Knudsen number) guarantees that the hydrodynamical 
description for a self-gravitating fluid becomes exact in the N, L — > oo limit for all scales ranging from 
the order L° till the order L. 

It must me noticed that the time between two collisions t co i = l/v ~ I is different both from 
the relaxation time and from the crossing time used in the literature. In particular, it is well known 

that EH ES| 

^crossinq 8 n AT 

^relaxation Jv 

This formula does not concern the time f co ; between two successive collisions. The time t co j is indeed 
very short due to the small angle behaviour of the gravitational cross section. For constant cross sections 
one finds a very different result for t co i [see ref. |16|]. 

In refs.EE] we thoroughly investigate the physics of the self-gravitating gas in thermal equilibrium. 
It is natural to study now the hydrodynamics of the self-gravitating fluid using p(r) = T p(f) as local 
equation of state. A first work on this direction is ref. |24| . 
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